02_Revisions
Suni Godbole
2024-02-16
1 Analysis File for manuscript
Read in Analytic Data files
pt.analytic.df <- read.csv(here::here("data", "PupilLightResponse_AnalyticSample_Demog.csv"),
header = T, stringsAsFactors = F)
pt.analytic.df$BMI_c <- pt.analytic.df$BMI - round(mean(pt.analytic.df$BMI), 2)
right_0.post.w <- read.csv(here::here("data", "PupilLightResponse_RightEye_Post_Wide.csv"),
header = T, stringsAsFactors = F)
rownames(right_0.post.w) <- right_0.post.w$subject_id
right_0.post <- read.csv(here::here("data", "PupilLightResponse_RightEye_Post_Long.csv"),
header = T, stringsAsFactors = F)2 Figure Individual Pupil Trajectories
right_0.post$user_cat2 <- ifelse(right_0.post$user_cat == "daily", 0,
ifelse(right_0.post$user_cat == "occasional", 1, 2))
right_0.post$user_cat2 <- factor(right_0.post$user_cat2,
levels = c(0,1,2),
labels = c("Daily Use", "Occasional Use", "No Use"))
post.user <- ggplot(right_0.post, aes(x=seconds, y = percent_change, group = subject_id))+
geom_line(show.legend = FALSE, alpha = 0.2)+
facet_grid(rows = vars(user_cat2))+
labs(y = "Percent change in pupil response",
x = "Seconds from start of light test")+
theme_bw()
post.user
2.0.1 Post use trajectors by use groups (individual plots) [manuscript]
jpeg(filename = here::here("figs", "Post_use_Indiv_Trajectories.jpeg"),
width = 6, height = 6, res = 300, units = "in")
post.user
dev.off()3 Analysis
3.1 SoFR for prediction of smoking
Another method that can be used to predict smoking status is a scalar-on-function regression model which uses differences between the whole trajectories to determine smoking status. However, when there are missing values (due to different test durations), all trajectories should be truncated to the shortest or the values should be imputed with fPCA. (See exploration of using SoFR line starting at 3000)
pctChg_vars <- names(right_0.post.w)[grepl("percent_change", names(right_0.post.w))]
pctChg_vars <- pctChg_vars[1:226]
sofr_impute_df <- right_0.post.w[, pctChg_vars]
rownames(sofr_impute_df) <- rownames(right_0.post.w)
sofr_fpca <- fpca.face(as.matrix(sofr_impute_df))
sofr_fpca2 <- sofr_fpca$Yhat
names(sofr_fpca2) <- names(sofr_impute_df)
## Adding in fpca imputed data where missing
sofr_imputed <- sofr_impute_df
sofr_imputed[is.na(sofr_imputed)] <- sofr_fpca2[is.na(sofr_imputed)]
ncols = ncol(sofr_imputed)
pt.analytic.df$female <- ifelse(pt.analytic.df$demo_gender == "Female", 1, 0)
smk.df <- pt.analytic.df[, c("subject_id", "user_cat", "female", "BMI_c")]
smk.df$smoker <- ifelse(pt.analytic.df$user_cat %in% c("daily", "occasional"), 1, 0)
sind <- seq(0, 1, len = ncols)
smat <- matrix(sind, nrow(smk.df), ncols, byrow = T)
# smk.df$Zsm <- I(Zsm)
smk.df$Zraw <- I(as.matrix(sofr_imputed))
smk.df$smat <- I(smat)
smk.df$lmat <- I(matrix(1/ncols, nrow(smk.df), ncols))
smk.df$zlmat <- I(smk.df$lmat * smk.df$Zraw)
# smk.df$bmiMat <- I(replicate(ncol(smk.df$smat), smk.df$BMI))
# smk.df$femaleMat <- I(replicate(ncol(smk.df$smat), smk.df$female))
smk_fglm_ps2 <- gam(smoker ~ s(smat, by=zlmat, bs = "cr", k = 30), data=smk.df,
method = "REML", family = binomial)
summary(smk_fglm_ps2)##
## Family: binomial
## Link function: logit
##
## Formula:
## smoker ~ s(smat, by = zlmat, bs = "cr", k = 30)
##
## Parametric coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.3826 0.5753 2.403 0.0162 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df Chi.sq p-value
## s(smat):zlmat 3.195 3.636 11.94 0.0189 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.141 Deviance explained = 13.7%
## -REML = 48.799 Scale est. = 1 n = 84pt.sofr.roc <- roc(smk_fglm_ps2$model$smoker, smk_fglm_ps2$fitted.values, ci = TRUE)
pt.sofr.roc.auc <- auc(pt.sofr.roc)
pt.sofr.auc.lci <- as.numeric(pt.sofr.roc$ci)[1]
pt.sofr.auc.uci <- as.numeric(pt.sofr.roc$ci)[3]3.2 Concurvity of SoFR model
concurvity(smk_fglm_ps2)## para s(smat):zlmat
## worst 0.8520874 0.8520874
## observed 0.8520874 0.3376077
## estimate 0.8520874 0.7471882concurvity(smk_fglm_ps2, full= FALSE)## $worst
## para s(smat):zlmat
## para 1.0000000 0.8520874
## s(smat):zlmat 0.8520874 1.0000000
##
## $observed
## para s(smat):zlmat
## para 1.0000000 0.3376077
## s(smat):zlmat 0.8520874 1.0000000
##
## $estimate
## para s(smat):zlmat
## para 1.0000000 0.7471882
## s(smat):zlmat 0.8520874 1.00000003.3 New model with sex and bmi as covariates
# Sex and gender added as scalar features
smk_fglm_cv <- gam(smoker ~ female + BMI_c +
s(smat, by=zlmat, bs = "cr", k = 30),
data=smk.df,
method = "REML", family = binomial)
summary(smk_fglm_cv)##
## Family: binomial
## Link function: logit
##
## Formula:
## smoker ~ female + BMI_c + s(smat, by = zlmat, bs = "cr", k = 30)
##
## Parametric coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.00598 0.69602 2.882 0.00395 **
## female -0.87859 0.52508 -1.673 0.09428 .
## BMI_c 0.06205 0.05924 1.047 0.29488
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df Chi.sq p-value
## s(smat):zlmat 3.065 3.464 12.02 0.0139 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.168 Deviance explained = 17.2%
## -REML = 48.252 Scale est. = 1 n = 84anova(smk_fglm_cv, smk_fglm_ps2, test = "Chisq")## Analysis of Deviance Table
##
## Model 1: smoker ~ female + BMI_c + s(smat, by = zlmat, bs = "cr", k = 30)
## Model 2: smoker ~ s(smat, by = zlmat, bs = "cr", k = 30)
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 77.138 89.662
## 2 78.922 93.456 -1.7843 -3.7937 0.12513.3.1 Extract Significant Regions in SoFR analysis
sigRegion <- function(df, crit.val, imp.var){
time <- imp.var[1]
est <- imp.var[2]
lci <- imp.var[3]
uci <- imp.var[4]
crit.rows <- which((df[, lci] < crit.val & df[, uci] < crit.val) |
(df[, lci] > crit.val & df[, uci] > crit.val))
# print(crit.rows)
group <- vector(mode = "numeric", length = length(crit.rows))
grp.ind <- 1
for(i in 1:length(crit.rows)){
if(i == 1){
group[i] <- grp.ind}
else(if(i > 1){
if(crit.rows[i] - crit.rows[i-1] == 1){
group[i] <- grp.ind
}
else(if(crit.rows[i] - crit.rows[i-1] != 1){
grp.ind <- grp.ind+1
group[i] <- grp.ind
})
})
}
# print(group)
num.grp <- unique(group)
ind.df <- data.frame(group = NA,
start = NA,
end = NA)
for(i in num.grp){
ind.df[i,] <- c(i, min(crit.rows[group == i]), max(crit.rows[group == i]))
}
# print(ind.df)
res.df <- data.frame(int.start = NA,
int.end = NA,
est.time = NA,
est = NA,
lci = NA,
uci = NA)
# print(nrow(ind.df))
for(i in 1:nrow(ind.df)){
sub.df <- df[(ind.df$start[i]):(ind.df$end[i]), ]
# print(sub.df[1, ])
if(sub.df[1, lci] > crit.val & sub.df[1, uci] >crit.val){
start1 <- round(sub.df[1, time], 2)
end1 <- round(sub.df[nrow(sub.df), time], 2)
or1 <- max(sub.df[ , est])
peaktime1 <- round(sub.df[sub.df[, est] == or1, time], 2)
lci1 <- round(sub.df[sub.df[, est] == or1, lci], 2)
uci1 <- round(sub.df[sub.df[, est] == or1, uci], 2)
or1 <- round(or1, 2)
res.df[i, ] <- c(start1, end1, peaktime1, or1, lci1, uci1)
# print(res.df)
}
if(sub.df[1, lci] < crit.val & sub.df[1, uci] < crit.val){
start2 <- round(sub.df[1, time], 2)
end2 <- round(sub.df[nrow(sub.df), time], 2)
or2 <- min(sub.df[, est])
peaktime2 <- round(sub.df[sub.df[, est] == or2, time], 2)
lci2 <- round(sub.df[sub.df[, est] == or2, lci], 2)
uci2 <- round(sub.df[sub.df[, est] == or2, uci], 2)
or2 <- round(or2, 2)
res.df[i, ] <- c(start2, end2, peaktime2, or2, lci2, uci2)
# print(res.df)
}
}
# crit.ind <- by(crit.rows, INDICES = list(group), FUN = function(x) c(min(x), max(x)))
return(res.df)
}## Covar Adjusted Model
df_sofr_sex <- data.frame("smat" = seq(0, 1, length.out = ncol(smk.df$smat)),
"zlmat" = 1, female = 0, BMI_c = 0)
sofr_pred_sex<- predict(smk_fglm_cv, newdata = df_sofr_sex, se.fit = TRUE, type = "iterms")
# Original Model (unadjusted)
df_sofr <- data.frame("smat" = seq(0, 1, length.out = ncol(smk.df$smat)),
"zlmat" = 1)
sofr_pred<- predict(smk_fglm_ps2, newdata = df_sofr, se.fit = TRUE, type = "iterms")
plot.sofr_sex <- data.frame("frame" = 0:(ncol(smk.df$smat) - 1),
"model_demog" = sofr_pred_sex$fit[,3],
"model_demog_se" = sofr_pred_sex$se.fit[,3],
"model_og" = sofr_pred$fit[,1],
"model_og_se" = sofr_pred$se.fit[,1])
plot.sofr_sex$lci <- exp(plot.sofr_sex$model_demog - 1.96*plot.sofr_sex$model_demog_se)
plot.sofr_sex$uci <- exp(plot.sofr_sex$model_demog + 1.96*plot.sofr_sex$model_demog_se)
plot.sofr_sex$OR <- exp(plot.sofr_sex$model_demog)
plot.sofr_sex$sec <- plot.sofr_sex$frame/30
plot.sofr_sex$p <- exp(plot.sofr_sex$model_demog)/(1+exp(plot.sofr_sex$model_demog))
sofr.sigRegions <- sigRegion(plot.sofr_sex, 1, imp.var = c("sec", "OR", "lci", "uci"))
sofr_plot_sex <- ggplot()+
geom_rect(data = sofr.sigRegions, aes(xmin = int.start, xmax = int.end,
ymin = -Inf, ymax = Inf), fill = "lightpink", alpha = 0.3)+
geom_hline(yintercept = 1, col = "darkorange")+
geom_line(data = plot.sofr_sex, aes(x=sec, y = exp(model_og)))+
geom_line(data = plot.sofr_sex, aes(x=sec, y = exp(model_og-1.96*model_og_se)), linetype = 'longdash')+
geom_line(data = plot.sofr_sex, aes(x=sec, y = exp(model_og+1.96*model_og_se)), linetype = 'longdash')+
geom_line(data = plot.sofr_sex, aes(x = sec, y = exp(model_demog)), col = 'purple')+
geom_line(data = plot.sofr_sex, aes(x=sec, y = exp(model_demog-1.96*model_demog_se)), linetype = 'longdash', col = "purple")+
geom_line(data = plot.sofr_sex, aes(x=sec, y = exp(model_demog+1.96*model_demog_se)), linetype = 'longdash', col = "purple")+
scale_x_continuous(expand = c(0, 0), limits = c(0,7.5),
breaks = c(0,2,4,5,6,7))+
scale_y_continuous(limits = c(0,5),
breaks = c(0,1,2,3,4,5)) +
labs(title = "",
x = "Seconds after start of light test",
y = "Odds ratio between smokers and non-smokers")+
theme_bw()
sofr_plot_sex
4 FoSR adjustment for covariates
4.1 Original model unadjusted for sex and BMI covariates.
right_0.gam.df <- merge(right_0.post, pt.analytic.df,
by = c("subject_id", "user_cat"))
right_0.gam.df$sind <- right_0.gam.df$frame/400
m.post_gam <- mgcv::gam(percent_change ~ s(sind, k=30, bs="cr") +
s(sind, by=occasional, k=30, bs = "cr") +
s(sind, by=daily, k=30, bs = "cr"),
data = right_0.gam.df, method = "REML")
## Create a matrix of residuals
post_gam.resid <- cbind(right_0.gam.df[!(is.na(right_0.gam.df$percent_change)),
c("subject_id", "frame")],
m.post_gam$residuals)
names(post_gam.resid) <- c("subject_id", "frame", "resid")
resid_mat <- reshape(post_gam.resid[, c("subject_id", "frame", "resid")],
timevar = "frame",
idvar = "subject_id",
direction = "wide")
rownames(resid_mat) <- resid_mat$subject_id
resid_mat$subject_id <- NULL
resid_mat <- as.matrix(resid_mat)
post_gam.fpca <- fpca.face(resid_mat, knots = 15)
#Create data frame of eigen functions and attach to original data
Phi_mat <- as.data.frame(post_gam.fpca$efunctions)
colnames(Phi_mat) <- paste0("Phi", seq(1, post_gam.fpca$npc))
Phi_mat$frame <- as.numeric(rownames(Phi_mat)) -1
right_0.gam.df <- merge(right_0.gam.df, Phi_mat,
by = "frame",
all.x = T)
right_0.gam.df <- right_0.gam.df[order(right_0.gam.df$subject_id, right_0.gam.df$frame), ]
right_0.gam.df$subject_id <- as.factor(right_0.gam.df$subject_id)
post_gam.fri <- mgcv::bam(percent_change ~
s(sind, k=30, bs="cr") +
s(sind, by=occasional, k=30, bs = "cr") +
s(sind, by=daily, k=30, bs = "cr") +
s(subject_id, by = Phi1, bs="re") + s(subject_id, by = Phi2, bs="re")+
s(subject_id, by = Phi3, bs="re") + s(subject_id, by = Phi4, bs="re"),
method = "fREML", data = right_0.gam.df, discrete = TRUE)
summary(post_gam.fri)##
## Family: gaussian
## Link function: identity
##
## Formula:
## percent_change ~ s(sind, k = 30, bs = "cr") + s(sind, by = occasional,
## k = 30, bs = "cr") + s(sind, by = daily, k = 30, bs = "cr") +
## s(subject_id, by = Phi1, bs = "re") + s(subject_id, by = Phi2,
## bs = "re") + s(subject_id, by = Phi3, bs = "re") + s(subject_id,
## by = Phi4, bs = "re")
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -10.3125 0.9197 -11.21 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(sind) 27.54 28.43 115.67 <2e-16 ***
## s(sind):occasional 18.63 21.95 70.14 <2e-16 ***
## s(sind):daily 18.62 21.95 35.42 <2e-16 ***
## s(subject_id):Phi1 82.23 84.00 23472.92 <2e-16 ***
## s(subject_id):Phi2 81.68 84.00 2775.02 <2e-16 ***
## s(subject_id):Phi3 82.23 84.00 887.87 <2e-16 ***
## s(subject_id):Phi4 81.05 84.00 1139.47 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.966 Deviance explained = 96.6%
## fREML = 63334 Scale est. = 2.6314 n = 326424.1.1 Add sex and BMI as covariates to the model (scalar and smoothed terms)
#original model: post_gam.fri
fosr_sex_gen.fri <- mgcv::bam(percent_change ~
female + BMI_c +
s(sind, k=30, bs="cr") +
s(sind, by=occasional, k=30, bs = "cr") +
s(sind, by=daily, k=30, bs = "cr") +
s(subject_id, by = Phi1, bs="re") + s(subject_id, by = Phi2, bs="re")+
s(subject_id, by = Phi3, bs="re") + s(subject_id, by = Phi4, bs="re"),
method = "fREML", data = right_0.gam.df, discrete = TRUE)
summary(fosr_sex_gen.fri)##
## Family: gaussian
## Link function: identity
##
## Formula:
## percent_change ~ female + BMI_c + s(sind, k = 30, bs = "cr") +
## s(sind, by = occasional, k = 30, bs = "cr") + s(sind, by = daily,
## k = 30, bs = "cr") + s(subject_id, by = Phi1, bs = "re") +
## s(subject_id, by = Phi2, bs = "re") + s(subject_id, by = Phi3,
## bs = "re") + s(subject_id, by = Phi4, bs = "re")
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -10.117542 0.922950 -10.962 <2e-16 ***
## female -0.206220 0.084028 -2.454 0.0141 *
## BMI_c -0.011684 0.009569 -1.221 0.2221
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(sind) 27.54 28.43 115.65 <2e-16 ***
## s(sind):occasional 18.63 21.95 70.40 <2e-16 ***
## s(sind):daily 18.62 21.95 35.37 <2e-16 ***
## s(subject_id):Phi1 82.23 84.00 23642.01 <2e-16 ***
## s(subject_id):Phi2 81.68 84.00 2758.50 <2e-16 ***
## s(subject_id):Phi3 82.23 84.00 884.39 <2e-16 ***
## s(subject_id):Phi4 81.05 84.00 1120.09 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.966 Deviance explained = 96.6%
## fREML = 63336 Scale est. = 2.6309 n = 32642anova(fosr_sex_gen.fri, post_gam.fri, test = "Chisq")## Analysis of Deviance Table
##
## Model 1: percent_change ~ female + BMI_c + s(sind, k = 30, bs = "cr") +
## s(sind, by = occasional, k = 30, bs = "cr") + s(sind, by = daily,
## k = 30, bs = "cr") + s(subject_id, by = Phi1, bs = "re") +
## s(subject_id, by = Phi2, bs = "re") + s(subject_id, by = Phi3,
## bs = "re") + s(subject_id, by = Phi4, bs = "re")
## Model 2: percent_change ~ s(sind, k = 30, bs = "cr") + s(sind, by = occasional,
## k = 30, bs = "cr") + s(sind, by = daily, k = 30, bs = "cr") +
## s(subject_id, by = Phi1, bs = "re") + s(subject_id, by = Phi2,
## bs = "re") + s(subject_id, by = Phi3, bs = "re") + s(subject_id,
## by = Phi4, bs = "re")
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 32238 84839
## 2 32240 84859 -2.003 -20.631 0.01989 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1df_pred_non_d <- data.frame("sind" = unique(right_0.gam.df$sind), "occasional" = 0, "daily" = 0,
"Phi1" = 0, "Phi2" = 0, "Phi3" = 0,
"Phi4" = 0, "female" = 0, "BMI_c" = 0,
"subject_id" = right_0.gam.df$subject_id[1])
df_pred_non <- data.frame("sind" = unique(right_0.gam.df$sind), "occasional" = 0, "daily" = 0,
"Phi1" = 0, "Phi2" = 0, "Phi3" = 0,
"Phi4" = 0,
"subject_id" = right_0.gam.df$subject_id[1])
pred_f1 <- predict(fosr_sex_gen.fri, newdata=df_pred_non_d, se.fit=TRUE, type = "response")
pred_f2 <- predict(post_gam.fri, newdata=df_pred_non, se.fit=TRUE, type = "response")
df_pred_occ_d <- data.frame("sind" = unique(right_0.gam.df$sind), "occasional" = 1, "daily" = 0,
"Phi1" = 0, "Phi2" = 0, "Phi3" = 0,
"Phi4" = 0, "female" = 0, "BMI_c" = 0,
"subject_id" = right_0.gam.df$subject_id[1])
df_pred_occ <- data.frame("sind" = unique(right_0.gam.df$sind), "occasional" = 1, "daily" = 0,
"Phi1" = 0, "Phi2" = 0, "Phi3" = 0,
"Phi4" = 0,
"subject_id" = right_0.gam.df$subject_id[1])
pred_f1_1 <- predict(fosr_sex_gen.fri, newdata=df_pred_occ_d, se.fit=TRUE, type = "response")
pred_f2_1 <- predict(post_gam.fri, newdata=df_pred_occ, se.fit=TRUE, type = "response")
df_pred_dly_d <- data.frame("sind" = unique(right_0.gam.df$sind), "occasional" = 0, "daily" = 1,
"Phi1" = 0, "Phi2" = 0, "Phi3" = 0,
"Phi4" = 0, "female" = 0, "BMI_c" = 0,
"subject_id" = right_0.gam.df$subject_id[1])
df_pred_dly <- data.frame("sind" = unique(right_0.gam.df$sind), "occasional" = 0, "daily" = 1,
"Phi1" = 0, "Phi2" = 0, "Phi3" = 0,
"Phi4" = 0,
"subject_id" = right_0.gam.df$subject_id[1])
pred_f1_2 <- predict(fosr_sex_gen.fri, newdata=df_pred_dly_d, se.fit=TRUE, type = "response")
pred_f2_2 <- predict(post_gam.fri, newdata=df_pred_dly, se.fit=TRUE, type = "response")
df <- data.frame(sec = 0:400/30,
user_group = c(rep("non-user", 802),
rep("occasional", 802),
rep("daily", 802)),
includ_demog = c(rep(1, 401), rep(0, 401),
rep(1, 401), rep(0, 401),
rep(1, 401), rep(0, 401)),
pct_chg = c(pred_f1$fit, pred_f2$fit, pred_f1_1$fit, pred_f2_1$fit, pred_f1_2$fit, pred_f2_2$fit))
df$includ_demog <- factor(df$includ_demog,
levels = c(0,1),
labels = c("Original Model",
"Model - Incl Demographics"))
demog_profiles <- ggplot(data = df, aes(x = sec, y = pct_chg, group = paste0(user_group, "_", includ_demog), col = user_group))+
geom_line(aes(linetype = as.factor(includ_demog)))+
labs(x = "Seconds from start of light test",
y = "Percent change in pupil response",
linetype = "Model Type",
col = "Use Group")Coefficient plots of demographic model vs original model
df_pred_occ <- data.frame("sind" = unique(right_0.gam.df$sind), "occasional" = 1, "daily" = 0,
"Phi1" = 0, "Phi2" = 0, "Phi3" = 0,
"Phi4" = 0, "female" = 0, "BMI_c" = 0,
"subject_id" = right_0.gam.df$subject_id[1])
df_pred_dly <- data.frame("sind" = unique(right_0.gam.df$sind), "occasional" = 0, "daily" = 1,
"Phi1" = 0, "Phi2" = 0, "Phi3" = 0,
"Phi4" = 0, "female" = 0, "BMI_c" = 0,
"subject_id" = right_0.gam.df$subject_id[1])
pred_f1 <- predict(fosr_sex_gen.fri, newdata=df_pred_occ, se.fit=TRUE, type = "terms")
pred_f2 <- predict(fosr_sex_gen.fri, newdata=df_pred_dly, se.fit=TRUE, type = "terms")
pred_sex_coef_df <- data.frame(seconds = seq(0, 400)/30,
occ_hat = pred_f1$fit[, 4],
occ_se = pred_f1$se.fit[, 4],
dly_hat = pred_f2$fit[, 5],
dly_se = pred_f2$se.fit[, 5])
pred_sex_coef_df$occ_lci <- pred_sex_coef_df$occ_hat - 2*pred_sex_coef_df$occ_se
pred_sex_coef_df$occ_uci <- pred_sex_coef_df$occ_hat + 2*pred_sex_coef_df$occ_se
pred_sex_coef_df$dly_lci <- pred_sex_coef_df$dly_hat - 2*pred_sex_coef_df$dly_se
pred_sex_coef_df$dly_uci <- pred_sex_coef_df$dly_hat + 2*pred_sex_coef_df$dly_se
occ_non.sigRegion <- sigRegion(pred_sex_coef_df, crit.val = 0, imp.var = c("seconds", "occ_hat", "occ_lci", "occ_uci"))
dly_non.sigRegion <- sigRegion(pred_sex_coef_df, crit.val = 0, imp.var = c("seconds", "dly_hat", "dly_lci", "dly_uci"))
post_occ_plt <- ggplot(data = pred_sex_coef_df, aes(x=seconds, y = occ_hat))+
geom_line()+
geom_line(aes(x = seconds, y = occ_hat + 2*occ_se), linetype = "longdash")+
geom_line(aes(x = seconds, y = occ_hat - 2*occ_se), linetype = "longdash")+
# geom_line(aes(x = seconds, y = 0), col = "darkred")+
geom_segment(data=occ_non.sigRegion, aes(x = int.start, y=0,
xend= int.end, yend = 0),
color = "darkred", linewidth =1.2)+
labs(title = "Occasional use vs No use",
y= "",
x = "")+
theme_bw()+
scale_x_continuous(expand = c(0, 0), limits = c(0,10),
breaks = c(0,2,4,6, 8, 10))+
scale_y_continuous(limits = c(-7, 7), breaks = c(-6,-4,-2,0,2,4,6))+
theme(rect = element_rect(fill = "transparent"))
post_dly_plt <- ggplot(data = pred_sex_coef_df, aes(x=seconds, y = dly_hat))+
geom_line()+
geom_line(aes(x = seconds, y = dly_hat + 2*dly_se), linetype = "longdash")+
geom_line(aes(x = seconds, y = dly_hat - 2*dly_se), linetype = "longdash")+
# geom_line(aes(x = seconds, y = 0), col = "darkred")+
geom_segment(data=dly_non.sigRegion, aes(x = int.start, y=0,
xend= int.end, yend = 0),
color = "darkred", linewidth =1.2)+
labs(title = "Daily use vs No use",
y = "",
x = "")+
theme_bw()+
scale_x_continuous(expand = c(0, 0), limits = c(0,10),
breaks = c(0,2,4,6, 8, 10))+
scale_y_continuous(limits = c(-7, 7), breaks = c(-6,-4,-2,0,2,4,6))+
theme(rect = element_rect(fill = "transparent"))
ylab <- textGrob(label = "Difference in percent change in pupil response", rot = 90,
gp = gpar(fontsize = 12))
posval <- textGrob(label = sapply(strwrap("Less Constriction in User", width = 15, simplify = F),
paste, collapse = "\n"), rot = 0,
gp = gpar(fontsize = 8))
negval <- textGrob(label = sapply(strwrap("More Constriction in User", width = 15, simplify = F),
paste, collapse = "\n"), rot = 0,
gp = gpar(fontsize = 8))
post_occ_coef <- grid.arrange(ylab, posval, negval, post_occ_plt,
ncol = 2, nrow = 2,
layout_matrix = rbind(c(1, 2, rep(4, 15)), c(1, 3, rep(4, 15))))
post_dly_coef <- grid.arrange(ylab, posval, negval, post_dly_plt,
ncol = 2, nrow = 2,
layout_matrix = rbind(c(1, 2, rep(4, 15)), c(1, 3, rep(4, 15))))
4.2 Review Comment Response Figure
demog.fosr <- grid.arrange(demog_profiles, post_occ_coef,
post_dly_coef,
layout_matrix = rbind(c(1,2),
c(3,4)))
demog.fosr## TableGrob (2 x 2) "arrange": 3 grobs
## z cells name grob
## 1 1 (1-1,1-1) arrange gtable[layout]
## 2 2 (1-1,2-2) arrange gtable[arrange]
## 3 3 (2-2,1-1) arrange gtable[arrange]